{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message:\n",
      "“\u001b[1m\u001b[22mThe `size` argument of `element_line()` is deprecated as of ggplot2 3.4.0.\n",
      "\u001b[36mℹ\u001b[39m Please use the `linewidth` argument instead.”\n",
      "Warning message:\n",
      "“\u001b[1m\u001b[22mThe `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.\n",
      "\u001b[36mℹ\u001b[39m Please use the `linewidth` argument instead.”\n"
     ]
    }
   ],
   "source": [
    "setwd(\"/frazer01/projects/PPC/analysis/ppc_eqtls\")\n",
    "\n",
    "source(\"scripts/packages.R\"  )\n",
    "source(\"scripts/input_data.R\")\n",
    "source(\"scripts/functions.R\" )\n",
    "source(\"scripts/coloc_functions.R\")\n",
    "source(\"scripts/3.4.coloc_adult/functions.R\")\n",
    "\n",
    "suppressMessages(library(coloc))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Objective: Check numbers in manuscript"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Abstract"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Total iPSC-PPC-unique: 1805\n",
      "\n"
     ]
    }
   ],
   "source": [
    "table = fread(\"reviews/tables/Table_S11_eQTL_Annotations.txt\", data.table = F)\n",
    "a = table %>% filter(category_annotation %in% c(\"ipsc_ppc-unique\", \"ipsc_ppc singleton\"))\n",
    "message(paste(\"Total iPSC-PPC-unique:\", nrow(a)))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "           \n",
       "            different partial same zero\n",
       "  different       109      68   45   46\n",
       "  partial          57     259  125   30\n",
       "  same             28     114    0    0\n",
       "  zero             76      86    0    0"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "\n",
       "alternative_splicing      gene_expression \n",
       "                 386                  657 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Total shared with adult: 1043 12.9 8108\n",
      "\n"
     ]
    }
   ],
   "source": [
    "a = table %>% filter(tissue == \"ipsc-ppc\" & category_annotation %in% c(\"fetal-adult\", \"fetal-islet\", \"fetal-pancreas\") &\n",
    "                            (islet_egene_overlap %in% c(\"partial\", \"different\") | pancreas_egene_overlap %in% c(\"partial\", \"different\")))  \n",
    "\n",
    "table(a$islet_egene_overlap, a$pancreas_egene_overlap)\n",
    "\n",
    "table(a$eqtl_phenotype)\n",
    "\n",
    "message(paste(\"Total shared with adult:\", nrow(a), \n",
    "              signif(nrow(a) / nrow(table[table$tissue == \"ipsc-ppc\",]) * 100, 3), \n",
    "              nrow(table[table$tissue == \"ipsc-ppc\",])))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Identification and characterization of gene and isoform eQTLs in fetal-like iPSC-PPCs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "gene_expression   isoform_usage \n",
       "           4065            4016 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "\n",
       "gene_expression   isoform_usage \n",
       "           4433            4232 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "5619\n",
      "\n",
      "Using 'Freq' as value column. Use 'value.var' to override\n",
      "\n",
      "Both: 1008\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "17.9391350774159"
      ],
      "text/latex": [
       "17.9391350774159"
      ],
      "text/markdown": [
       "17.9391350774159"
      ],
      "text/plain": [
       "[1] 17.93914"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Only gene exp: 3057\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "54.4046983449012"
      ],
      "text/latex": [
       "54.4046983449012"
      ],
      "text/markdown": [
       "54.4046983449012"
      ],
      "text/plain": [
       "[1] 54.4047"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Only isof use: 1554\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "27.6561665776829"
      ],
      "text/latex": [
       "27.6561665776829"
      ],
      "text/markdown": [
       "27.6561665776829"
      ],
      "text/plain": [
       "[1] 27.65617"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "table = fread(\"reviews/tables/Table_S7_LeadSNP.txt\", data.table = F) %>% distinct()\n",
    "\n",
    "# Number of eGenes and eIsoforms\n",
    "a = table %>% filter(discovery_order == 0 & egene == T)\n",
    "table(a$eqtl_phenotype)\n",
    "\n",
    "# Number of egQTls and eiQTLs\n",
    "a = table %>% filter(egene == T)\n",
    "table(a$eqtl_phenotype)\n",
    "\n",
    "# Number of genes affected by either gene expression or isoform usage\n",
    "a = table %>% filter(egene == T)\n",
    "message(length(unique(a$gene_id)))\n",
    "all = length(unique(a$gene_id))\n",
    "\n",
    "b = suppressWarnings(dcast(data.frame(table(a$gene_id, a$eqtl_phenotype)), Var1 ~ Var2))\n",
    "\n",
    "message(paste(\"Both:\", length(unique(b[b$gene_expression != 0 & b$isoform_usage != 0,]$Var1)) ))\n",
    "length(unique(b[b$gene_expression != 0 & b$isoform_usage != 0,]$Var1)) / all * 100\n",
    "\n",
    "message(paste(\"Only gene exp:\", length(unique(b[b$gene_expression != 0 & b$isoform_usage == 0,]$Var1)) ))\n",
    "length(unique(b[b$gene_expression != 0 & b$isoform_usage == 0,]$Var1)) / all * 100\n",
    "\n",
    "message(paste(\"Only isof use:\", length(unique(b[b$gene_expression == 0 & b$isoform_usage != 0,]$Var1)) ))\n",
    "length(unique(b[b$gene_expression == 0 & b$isoform_usage != 0,]$Var1)) / all * 100\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       " H3  H4 \n",
       "147 671 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. \n",
       " 0.8001  0.8998  0.9565  0.9373  0.9866  1.0000 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "# genes with H3 or H4 with their isoforms: 410\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "40.6746031746032"
      ],
      "text/latex": [
       "40.6746031746032"
      ],
      "text/markdown": [
       "40.6746031746032"
      ],
      "text/plain": [
       "[1] 40.6746"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "# genes with only H3: 333\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "81.219512195122"
      ],
      "text/latex": [
       "81.219512195122"
      ],
      "text/markdown": [
       "81.219512195122"
      ],
      "text/plain": [
       "[1] 81.21951"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "# genes with only H3: 38\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "9.26829268292683"
      ],
      "text/latex": [
       "9.26829268292683"
      ],
      "text/markdown": [
       "9.26829268292683"
      ],
      "text/plain": [
       "[1] 9.268293"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "# genes with both H3 and H4: 39\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "9.51219512195122"
      ],
      "text/latex": [
       "9.51219512195122"
      ],
      "text/markdown": [
       "9.51219512195122"
      ],
      "text/plain": [
       "[1] 9.512195"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "table = fread(\"reviews/tables/Table_S10_PPC_Gene_Isoform_Colocalization.txt\", data.table = F)\n",
    "\n",
    "a = table %>% filter(max_model_pp >= 0.8 & likely_model %in% c(\"H3\", \"H4\"))\n",
    "table(a$likely_model)\n",
    "summary(a$max_model_pp)\n",
    "\n",
    "message(paste(\"# genes with H3 or H4 with their isoforms:\", length(unique(a$gene_id.1)) ))\n",
    "length(unique(a$gene_id.1)) / 1008 * 100\n",
    "\n",
    "b = data.frame(table(a$gene_id.1, a$likely_model))\n",
    "b = suppressMessages(suppressWarnings(dcast(b, Var1 ~ Var2)))\n",
    "\n",
    "all = length(unique(a$gene_id.1))\n",
    "\n",
    "message(paste(\"# genes with only H3:\", length(unique(b[b$H3 == 0 & b$H4 != 0,]$Var1)) ))\n",
    "length(unique(b[b$H3 == 0 & b$H4 != 0,]$Var1)) / all * 100\n",
    "\n",
    "message(paste(\"# genes with only H3:\", length(unique(b[b$H3 != 0 & b$H4 == 0,]$Var1)) ))\n",
    "length(unique(b[b$H3 != 0 & b$H4 == 0,]$Var1)) / all * 100\n",
    "\n",
    "message(paste(\"# genes with both H3 and H4:\", length(unique(b[b$H3 != 0 & b$H4 != 0,]$Var1)) ))\n",
    "length(unique(b[b$H3 != 0 & b$H4 != 0,]$Var1)) / all * 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Most fetal-like and adult islet eGenes show developmental stage specificity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"dataframe\">\n",
       "<caption>A data.frame: 2 × 24</caption>\n",
       "<thead>\n",
       "\t<tr><th></th><th scope=col>eqtl_id.1</th><th scope=col>eqtl_id.2</th><th scope=col>transcript_id.1</th><th scope=col>transcript_id.2</th><th scope=col>gene_id.1</th><th scope=col>gene_id.2</th><th scope=col>gene_name.1</th><th scope=col>gene_name.2</th><th scope=col>eqtl_phenotype.1</th><th scope=col>eqtl_phenotype.2</th><th scope=col>⋯</th><th scope=col>nsnps</th><th scope=col>PP.H0.abf</th><th scope=col>PP.H1.abf</th><th scope=col>PP.H2.abf</th><th scope=col>PP.H3.abf</th><th scope=col>PP.H4.abf</th><th scope=col>likely_model</th><th scope=col>max_model_pp</th><th scope=col>topsnp</th><th scope=col>topsnp_pp</th></tr>\n",
       "\t<tr><th></th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>⋯</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;dbl&gt;</th></tr>\n",
       "</thead>\n",
       "<tbody>\n",
       "\t<tr><th scope=row>1</th><td>ipsc-ppc_0_ENSG00000003249</td><td>islet_0_ENSG00000003249</td><td>ENSG00000003249</td><td>ENSG00000003249</td><td>ENSG00000003249</td><td>ENSG00000003249</td><td>DBNDD1</td><td>DBNDD1</td><td>gene_expression</td><td>gene_expression</td><td>⋯</td><td>1847</td><td>6.764578e-05</td><td>4.952350e-05</td><td>0.07052776</td><td>0.05075480</td><td>0.8786003</td><td>H4</td><td>0.8786003</td><td>VAR_16_90081660_A_G</td><td>0.78119063</td></tr>\n",
       "\t<tr><th scope=row>2</th><td>ipsc-ppc_0_ENSG00000004534</td><td>islet_0_ENSG00000004534</td><td>ENSG00000004534</td><td>ENSG00000004534</td><td>ENSG00000004534</td><td>ENSG00000004534</td><td>RBM6  </td><td>RBM6  </td><td>gene_expression</td><td>gene_expression</td><td>⋯</td><td>1152</td><td>1.484825e-07</td><td>6.369398e-07</td><td>0.02309440</td><td>0.09818851</td><td>0.8787163</td><td>H4</td><td>0.8787163</td><td>VAR_3_50174848_G_A </td><td>0.07459922</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "A data.frame: 2 × 24\n",
       "\\begin{tabular}{r|lllllllllllllllllllll}\n",
       "  & eqtl\\_id.1 & eqtl\\_id.2 & transcript\\_id.1 & transcript\\_id.2 & gene\\_id.1 & gene\\_id.2 & gene\\_name.1 & gene\\_name.2 & eqtl\\_phenotype.1 & eqtl\\_phenotype.2 & ⋯ & nsnps & PP.H0.abf & PP.H1.abf & PP.H2.abf & PP.H3.abf & PP.H4.abf & likely\\_model & max\\_model\\_pp & topsnp & topsnp\\_pp\\\\\n",
       "  & <chr> & <chr> & <chr> & <chr> & <chr> & <chr> & <chr> & <chr> & <chr> & <chr> & ⋯ & <int> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <chr> & <dbl> & <chr> & <dbl>\\\\\n",
       "\\hline\n",
       "\t1 & ipsc-ppc\\_0\\_ENSG00000003249 & islet\\_0\\_ENSG00000003249 & ENSG00000003249 & ENSG00000003249 & ENSG00000003249 & ENSG00000003249 & DBNDD1 & DBNDD1 & gene\\_expression & gene\\_expression & ⋯ & 1847 & 6.764578e-05 & 4.952350e-05 & 0.07052776 & 0.05075480 & 0.8786003 & H4 & 0.8786003 & VAR\\_16\\_90081660\\_A\\_G & 0.78119063\\\\\n",
       "\t2 & ipsc-ppc\\_0\\_ENSG00000004534 & islet\\_0\\_ENSG00000004534 & ENSG00000004534 & ENSG00000004534 & ENSG00000004534 & ENSG00000004534 & RBM6   & RBM6   & gene\\_expression & gene\\_expression & ⋯ & 1152 & 1.484825e-07 & 6.369398e-07 & 0.02309440 & 0.09818851 & 0.8787163 & H4 & 0.8787163 & VAR\\_3\\_50174848\\_G\\_A  & 0.07459922\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "A data.frame: 2 × 24\n",
       "\n",
       "| <!--/--> | eqtl_id.1 &lt;chr&gt; | eqtl_id.2 &lt;chr&gt; | transcript_id.1 &lt;chr&gt; | transcript_id.2 &lt;chr&gt; | gene_id.1 &lt;chr&gt; | gene_id.2 &lt;chr&gt; | gene_name.1 &lt;chr&gt; | gene_name.2 &lt;chr&gt; | eqtl_phenotype.1 &lt;chr&gt; | eqtl_phenotype.2 &lt;chr&gt; | ⋯ ⋯ | nsnps &lt;int&gt; | PP.H0.abf &lt;dbl&gt; | PP.H1.abf &lt;dbl&gt; | PP.H2.abf &lt;dbl&gt; | PP.H3.abf &lt;dbl&gt; | PP.H4.abf &lt;dbl&gt; | likely_model &lt;chr&gt; | max_model_pp &lt;dbl&gt; | topsnp &lt;chr&gt; | topsnp_pp &lt;dbl&gt; |\n",
       "|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|\n",
       "| 1 | ipsc-ppc_0_ENSG00000003249 | islet_0_ENSG00000003249 | ENSG00000003249 | ENSG00000003249 | ENSG00000003249 | ENSG00000003249 | DBNDD1 | DBNDD1 | gene_expression | gene_expression | ⋯ | 1847 | 6.764578e-05 | 4.952350e-05 | 0.07052776 | 0.05075480 | 0.8786003 | H4 | 0.8786003 | VAR_16_90081660_A_G | 0.78119063 |\n",
       "| 2 | ipsc-ppc_0_ENSG00000004534 | islet_0_ENSG00000004534 | ENSG00000004534 | ENSG00000004534 | ENSG00000004534 | ENSG00000004534 | RBM6   | RBM6   | gene_expression | gene_expression | ⋯ | 1152 | 1.484825e-07 | 6.369398e-07 | 0.02309440 | 0.09818851 | 0.8787163 | H4 | 0.8787163 | VAR_3_50174848_G_A  | 0.07459922 |\n",
       "\n"
      ],
      "text/plain": [
       "  eqtl_id.1                  eqtl_id.2               transcript_id.1\n",
       "1 ipsc-ppc_0_ENSG00000003249 islet_0_ENSG00000003249 ENSG00000003249\n",
       "2 ipsc-ppc_0_ENSG00000004534 islet_0_ENSG00000004534 ENSG00000004534\n",
       "  transcript_id.2 gene_id.1       gene_id.2       gene_name.1 gene_name.2\n",
       "1 ENSG00000003249 ENSG00000003249 ENSG00000003249 DBNDD1      DBNDD1     \n",
       "2 ENSG00000004534 ENSG00000004534 ENSG00000004534 RBM6        RBM6       \n",
       "  eqtl_phenotype.1 eqtl_phenotype.2 ⋯ nsnps PP.H0.abf    PP.H1.abf   \n",
       "1 gene_expression  gene_expression  ⋯ 1847  6.764578e-05 4.952350e-05\n",
       "2 gene_expression  gene_expression  ⋯ 1152  1.484825e-07 6.369398e-07\n",
       "  PP.H2.abf  PP.H3.abf  PP.H4.abf likely_model max_model_pp topsnp             \n",
       "1 0.07052776 0.05075480 0.8786003 H4           0.8786003    VAR_16_90081660_A_G\n",
       "2 0.02309440 0.09818851 0.8787163 H4           0.8787163    VAR_3_50174848_G_A \n",
       "  topsnp_pp \n",
       "1 0.78119063\n",
       "2 0.07459922"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "H4: 701\n",
      "\n",
      "H3: 94\n",
      "\n"
     ]
    }
   ],
   "source": [
    "table = fread(\"reviews/tables/Table_S10_PPC_Islets_Colocalization_eGene.txt\", data.table = F)\n",
    "head(table,2)\n",
    "\n",
    "# Genes with at least one H4\n",
    "message(paste(\"H4:\", length(unique(table[table$likely_model == \"H4\",]$gene_id.1))))\n",
    "\n",
    "# Genes with only H3 and no H4\n",
    "a = table[table$likely_model == \"H4\",]\n",
    "message(paste(\"H3:\", length(unique(table[table$likely_model == \"H3\" & !table$gene_id.1 %in% a$gene_id.1,]$gene_id.1))))\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Developmental stage-unique and shared egQTLs "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "ipsc_ppc singleton    islet singleton pancreas singleton \n",
       "               887                703               1927 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "table = fread(\"reviews/tables/Table_S11_eQTL_Annotations.txt\", data.table = F) \n",
    "\n",
    "a = table %>% filter(eqtl_phenotype %like% \"gene\" & category_annotation %like% \"singleton\")\n",
    "table(a$category_annotation)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "'data.frame':\t2720 obs. of  12 variables:\n",
      " $ module_id             : chr  \"GE_1_1\" \"GE_1_10\" \"GE_1_100\" \"GE_1_101\" ...\n",
      " $ associations          : chr  \"ipsc-ppc_0_ENSG00000116874,ipsc-ppc_0_ENSG00000231365,islet_0_ENSG00000116874,pancreas_0_ENSG00000231365,pancre\"| __truncated__ \"islet_0_ENSG00000007341,islet_0_ENSG00000134245,pancreas_0_ENSG00000007341,pancreas_0_ENSG00000116489,pancreas_\"| __truncated__ \"pancreas_0_ENSG00000225217,pancreas_0_ENSG00000244682\" \"pancreas_0_ENSG00000226026,pancreas_0_ENSG00000228852,ipsc-ppc_0_ENSG00000235501\" ...\n",
      " $ number_assocs         : int  5 6 2 3 2 2 3 2 2 2 ...\n",
      " $ number_ipsc_ppc_assocs: int  2 1 0 1 0 0 0 0 0 0 ...\n",
      " $ number_islet_assocs   : int  1 2 0 0 0 0 0 0 0 0 ...\n",
      " $ number_pancreas_assocs: int  2 3 2 2 2 2 3 2 2 2 ...\n",
      " $ islet_egene_overlap   : chr  \"partial\" \"partial\" \"\" \"zero\" ...\n",
      " $ pancreas_egene_overlap: chr  \"same\" \"partial\" \"\" \"different\" ...\n",
      " $ egene_overlap_category: chr  \"C\" \"C\" \"\" \"D\" ...\n",
      " $ module_pass           : logi  FALSE TRUE TRUE TRUE TRUE TRUE ...\n",
      " $ category_annotation   : chr  \"module_failed\" \"fetal-adult\" \"pancreas-unique\" \"fetal-pancreas\" ...\n",
      " $ notes                 : chr  \"Failed module requirements\" \"N/A\" \"N/A\" \"N/A\" ...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "# modules passed criteria: 1852\n",
      "\n",
      "# mod, 2 egQTLs: 939 50.7\n",
      "\n",
      "# mod, > 2 egQTLs: 913 49.3\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. \n",
       "  3.000   3.000   3.000   3.954   4.000  20.000 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "table = fread(\"reviews/tables/Table_S12_ModuleSummary_Gene.txt\", data.table = F)\n",
    "table2 = fread(\"reviews/tables/Table_S11_eQTL_Annotations.txt\")\n",
    "str(table)\n",
    "\n",
    "a = table %>% filter(module_pass == T & !category_annotation %like% \"ambi\")\n",
    "message(paste(\"# modules passed criteria:\", length(unique(a$module_id)) ))\n",
    "\n",
    "b = a %>% filter(number_assocs == 2)\n",
    "message(paste(\"# mod, 2 egQTLs:\", length(unique(b$module_id)), signif(length(unique(b$module_id)) / length(unique(a$module_id)) * 100, 4)))\n",
    "\n",
    "b = a %>% filter(number_assocs > 2)\n",
    "message(paste(\"# mod, > 2 egQTLs:\", length(unique(b$module_id)), signif(length(unique(b$module_id)) / length(unique(a$module_id)) * 100, 4)))\n",
    "\n",
    "summary(a[a$number_assocs > 2,]$number_assocs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "# mod tissue-unique: 199 10.7 1852\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "\n",
       "   adult-shared     fetal-adult     fetal-islet  fetal-pancreas ipsc_ppc-unique \n",
       "            670             652              53             278              10 \n",
       "   islet-unique pancreas-unique \n",
       "             30             159 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "\n",
       "ipsc-ppc    islet pancreas \n",
       "      10       30      159 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "\n",
       "ipsc-ppc    islet pancreas \n",
       "      21       62      354 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "b = a[a$category_annotation %like% \"unique\",] \n",
    "message(paste(\"# mod tissue-unique:\", nrow(b), signif(nrow(b) / nrow(a) * 100, 3), nrow(a)))\n",
    "table(a$category_annotation)\n",
    "\n",
    "b = table2 %>% filter(category_annotation %like% \"unique\" & eqtl_phenotype == \"gene_expression\") %>% select(module_id, tissue) %>% distinct()\n",
    "table(b$tissue)\n",
    "\n",
    "b = table2 %>% filter(category_annotation %like% \"unique\" & eqtl_phenotype == \"gene_expression\") %>% distinct()\n",
    "table(b$tissue)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "# mod tissue-sharing: 1653 89.25 1852\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "\n",
       "  adult-shared    fetal-adult    fetal-islet fetal-pancreas \n",
       "           670            652             53            278 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "b = a[!a$category_annotation %like% \"unique\" & a$category_annotation != \"ambiguous\",] \n",
    "message(paste(\"# mod tissue-sharing:\", nrow(b), signif(nrow(b) / nrow(a) * 100, 4), nrow(a)))\n",
    "table(b$category_annotation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "  adult-shared    fetal-adult    fetal-islet fetal-pancreas \n",
       "           670            652             53            278 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "# total fetal-adult sharing: 983\n",
      "\n"
     ]
    }
   ],
   "source": [
    "a = table2 %>% \n",
    "    filter(!category_annotation %like% \"failed\" & \n",
    "           !category_annotation %like% \"singleton\" & \n",
    "           !category_annotation %like% \"unique\" & \n",
    "           category_annotation != \"ambiguous\" & \n",
    "           eqtl_phenotype == \"gene_expression\") %>% \n",
    "    select(module_id, category_annotation) %>% distinct()\n",
    "table(a$category_annotation)\n",
    "\n",
    "message(paste(\"# total fetal-adult sharing:\", nrow(a[a$category_annotation %like% \"fetal\",])))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "   fetal-adult    fetal-islet fetal-pancreas \n",
       "           652             53            278 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "# mod fetal-adult sharing: 983 36.14 2720\n",
      "\n",
      "how many eQTLs:\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "\n",
       "ipsc-ppc    islet pancreas \n",
       "    1122      870     1394 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "a = table[!table$category_annotation %like% \"failed\" & \n",
    "        table$category_annotation != \"adult-shared\" & \n",
    "        !table$category_annotation %like% \"unique\" & \n",
    "        table$category_annotation != \"ambiguous\",] \n",
    "table(a$category_annotation)\n",
    "\n",
    "message(paste(\"# mod fetal-adult sharing:\", nrow(a), signif(nrow(a) / nrow(table) * 100, 4), nrow(table)))\n",
    "\n",
    "a = table2 %>% filter(category_annotation %in% c(\"fetal-adult\", \"fetal-islets\", \"fetal-pancreas\") & eqtl_phenotype == \"gene_expression\") %>% \n",
    "    select(eqtl_id, tissue) %>% distinct()\n",
    "message(paste(\"how many eQTLs:\"))\n",
    "table(a$tissue)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### conclusion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "# total ppc unique: 1805\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "                      \n",
       "                       combinatorial singleton\n",
       "  alternative_splicing           266       631\n",
       "  gene_expression                 21       887"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "\n",
       "combinatorial     singleton \n",
       "          287          1518 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "a = table2 %>% filter(category_annotation %in% c(\"ipsc_ppc singleton\" , \"ipsc_ppc-unique\"))\n",
    "message(paste(\"# total ppc unique:\", length(unique(a$eqtl_id))))\n",
    "\n",
    "\n",
    "table(a$eqtl_phenotype, a$eqtl_type)\n",
    "\n",
    "table(a$eqtl_type)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "# total ppc shared with adult: 1977\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "\n",
       "alternative_splicing      gene_expression \n",
       "                 802                 1175 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# iPSC-PPC eQTLs shared with adult\n",
    "a = table2 %>% filter(category_annotation %in% c(\"fetal-adult\", \"fetal-islet\", \"fetal-pancreas\") & tissue == \"ipsc-ppc\")\n",
    "message(paste(\"# total ppc shared with adult:\", nrow(a)))\n",
    "table(a$eqtl_phenotype)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "# failed: 4326\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "\n",
       "alternative_splicing      gene_expression \n",
       "                2260                 2066 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "a = table2 %>% filter(category_annotation %in% c(\"module_failed\", \"ambiguous\") & tissue == \"ipsc-ppc\")\n",
    "message(paste(\"# failed:\", nrow(a)))\n",
    "table(a$eqtl_phenotype)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Regulatory plasticity in combinatorial egQTLs shared between fetal-like and adult pancreatic tissues"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 128,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "# total modules shared b/t ppc and adult: 983\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "\n",
       "   fetal-adult    fetal-islet fetal-pancreas \n",
       "           652             53            278 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "\n",
       "  A   B   C   D   E \n",
       "200 305 350  88  40 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "\n",
       "    A     B     C     D     E \n",
       "20.30 31.00 35.60  8.95  4.07 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "# same genes: 505\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "51.3733468972533"
      ],
      "text/latex": [
       "51.3733468972533"
      ],
      "text/markdown": [
       "51.3733468972533"
      ],
      "text/plain": [
       "[1] 51.37335"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "# same genes: 478\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "48.6266531027467"
      ],
      "text/latex": [
       "48.6266531027467"
      ],
      "text/markdown": [
       "48.6266531027467"
      ],
      "text/plain": [
       "[1] 48.62665"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "table = fread(\"reviews/tables/Table_S12_ModuleSummary_Gene.txt\", data.table = F) %>% distinct()\n",
    "table2 = fread(\"reviews/tables/Table_S11_eQTL_Annotations.txt\", data.table = F)\n",
    "\n",
    "a = table %>% filter(category_annotation %in% c(\"fetal-adult\", \"fetal-islet\", \"fetal-pancreas\"))\n",
    "message(paste(\"# total modules shared b/t ppc and adult:\", nrow(a)))\n",
    "table(a$category_annotation)\n",
    "table(a$egene_overlap_category)\n",
    "\n",
    "signif(table(a$egene_overlap_category) / nrow(a) * 100, 3)\n",
    "\n",
    "b = a %>% filter(egene_overlap_category %in% c(\"A\", \"B\"))\n",
    "message(paste(\"# same genes:\", nrow(b)))\n",
    "nrow(b) / nrow(a) * 100\n",
    "\n",
    "b = a %>% filter(egene_overlap_category %in% c(\"C\", \"D\", \"E\"))\n",
    "message(paste(\"# same genes:\", nrow(b)))\n",
    "nrow(b) / nrow(a) * 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Another way to get category modules (sanity-check)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 129,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "# mod in A: 200 20.35\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. \n",
       "   1.00    1.00    1.00    1.01    1.00    2.00 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "b = a %>% filter( (islet_egene_overlap == \"zero\" & pancreas_egene_overlap == \"same\") | (islet_egene_overlap == \"same\" & pancreas_egene_overlap == \"zero\") )\n",
    "message(paste(\"# mod in A:\", nrow(b), signif(nrow(b) / nrow(a) * 100, 4)))\n",
    "\n",
    "b = table2 %>% filter(eqtl_phenotype == \"gene_expression\" & module_id %in% b$module_id) %>% select(module_id, gene_id) %>% distinct()\n",
    "summary(data.frame(table(b$module_id))$Freq)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 131,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "# mod in B: 305 31.03\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. \n",
       "   1.00    1.00    1.00    1.01    1.00    2.00 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "b = a %>% filter( (islet_egene_overlap == \"same\" & pancreas_egene_overlap == \"same\") | \n",
    "                  (islet_egene_overlap == \"same\" & pancreas_egene_overlap == \"same\") )\n",
    "message(paste(\"# mod in B:\", nrow(b), signif(nrow(b) / nrow(a) * 100, 4)))\n",
    "\n",
    "b = table2 %>% filter(eqtl_phenotype == \"gene_expression\" & module_id %in% b$module_id) %>% select(module_id, gene_id) %>% distinct()\n",
    "summary(data.frame(table(b$module_id))$Freq)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "# mod in C: 350 35.61\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. \n",
       "  2.000   2.000   2.000   2.831   3.000  12.000 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "b = a %>% filter(islet_egene_overlap %like% \"partial\" | pancreas_egene_overlap %like% \"partial\")\n",
    "message(paste(\"# mod in C:\", nrow(b), signif(nrow(b) / nrow(a) * 100, 4)))\n",
    "\n",
    "b = table2 %>% filter(eqtl_phenotype == \"gene_expression\" & module_id %in% b$module_id) %>% select(module_id, gene_id) %>% distinct()\n",
    "summary(data.frame(table(b$module_id))$Freq)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 133,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "# mod in D: 88 9\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. \n",
       "  2.000   2.000   2.000   2.284   2.000   5.000 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "b = a %>% filter( (islet_egene_overlap == \"different\" & pancreas_egene_overlap %in% c(\"zero\", \"same\")) | \n",
    "                  (islet_egene_overlap %in% c(\"zero\", \"same\") & pancreas_egene_overlap == \"different\") )\n",
    "message(paste(\"# mod in D:\", nrow(b), signif(nrow(b) / nrow(a) * 100, 2)))\n",
    "\n",
    "b = table2 %>% filter(eqtl_phenotype == \"gene_expression\" & module_id %in% b$module_id) %>% select(module_id, gene_id) %>% distinct()\n",
    "summary(data.frame(table(b$module_id))$Freq)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 135,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "# mod in E: 40 4.1\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. \n",
       "   2.00    2.00    3.00    3.15    4.00    7.00 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "b = a %>% filter( islet_egene_overlap == \"different\" & pancreas_egene_overlap == \"different\" )\n",
    "message(paste(\"# mod in E:\", nrow(b), signif(nrow(b) / nrow(a) * 100, 2)))\n",
    "\n",
    "b = table2 %>% filter(eqtl_phenotype == \"gene_expression\" & module_id %in% b$module_id) %>% select(module_id, gene_id) %>% distinct()\n",
    "summary(data.frame(table(b$module_id))$Freq)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 140,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "total modules shared b/t ppc and adult: 655\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "\n",
       "ipsc-ppc    islet pancreas \n",
       "    1043      934     1111 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "mod1 = fread(\"reviews/tables/Table_S12_ModuleSummary_Gene.txt\", data.table = F) %>% filter(!category_annotation %like% \"module_fail\") %>% distinct()\n",
    "mod2 = fread(\"reviews/tables/Table_S12_ModuleSummary_AltSplicing.txt\", data.table = F) %>% filter(!category_annotation %like% \"module_fail\") %>% distinct()\n",
    "mod = rbind(mod1, mod2) %>% filter(category_annotation %in% c(\"fetal-adult\", \"fetal-islet\", \"fetal-pancreas\") & egene_overlap_category %in% c(\"C\", \"D\", \"E\")) \n",
    "\n",
    "message(paste(\"total modules shared b/t ppc and adult:\", nrow(mod)))\n",
    "\n",
    "table(table2[table2$module_id %in% mod$module_id,]$tissue)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Associations of developmental stage-unique eQTLs with pancreatic traits and disease phenotypes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 141,
   "metadata": {},
   "outputs": [],
   "source": [
    "gwas_coloc = fread(\"reviews/tables/Table_S14_GWAS_Colocalization.txt\", data.table = F)\n",
    "eqtl_annot = fread(\"reviews/tables/Table_S11_eQTL_Annotations.txt\", data.table = F)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Singleton"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 145,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "# total singletons: 6101\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "\n",
       "alternative_splicing      gene_expression \n",
       "                2584                 3517 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "a = eqtl_annot %>% filter(category_annotation %like% \"singleton\")\n",
    "message(paste(\"# total singletons:\", length(unique(a$eqtl_id))))\n",
    "table(a$eqtl_phenotype)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 148,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "# singletons with GWAS coloc: 118 1.934\n",
      "\n",
      "# ppc-singleton: 21 1518 1.383\n",
      "\n",
      "# islet-singleton: 57 2225 2.562\n",
      "\n",
      "# panc-singleton: 40 2358 1.696\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "\n",
       "ipsc-ppc    islet pancreas \n",
       "      21       57       40 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "b = gwas_coloc[gwas_coloc$eqtl_id %in% a$eqtl_id,]\n",
    "message(paste(\"# singletons with GWAS coloc:\", length(unique(b$eqtl_id)), signif(length(unique(b$eqtl_id)) / length(unique(a$eqtl_id)) * 100, 4) ))\n",
    "message(paste(\"# ppc-singleton:\", length(unique(b[b$tissue == \"ipsc-ppc\",]$eqtl_id)), length(unique(a[a$tissue == \"ipsc-ppc\",]$eqtl_id)), signif(length(unique(b[b$tissue == \"ipsc-ppc\",]$eqtl_id)) / length(unique(a[a$tissue == \"ipsc-ppc\",]$eqtl_id)) * 100, 4) ))\n",
    "message(paste(\"# islet-singleton:\", length(unique(b[b$tissue == \"islet\",]$eqtl_id)), length(unique(a[a$tissue == \"islet\",]$eqtl_id)), signif(length(unique(b[b$tissue == \"islet\",]$eqtl_id)) / length(unique(a[a$tissue == \"islet\",]$eqtl_id)) * 100, 4) ))\n",
    "message(paste(\"# panc-singleton:\", length(unique(b[b$tissue == \"pancreas\",]$eqtl_id)), length(unique(a[a$tissue == \"pancreas\",]$eqtl_id)), signif(length(unique(b[b$tissue == \"pancreas\",]$eqtl_id)) / length(unique(a[a$tissue == \"pancreas\",]$eqtl_id)) * 100, 4) ))\n",
    "\n",
    "table(unique(b[,c(\"eqtl_id\", \"tissue\")])$tissue)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 154,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "eQTls colocalized with > 1 traits: 38\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. \n",
       "  2.000   2.000   2.000   2.711   3.000   6.000 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "b = gwas_coloc[gwas_coloc$eqtl_id %in% a$eqtl_id,] %>% select(eqtl_id, trait_id) %>% distinct()\n",
    "b = data.frame(table(b$eqtl_id)) %>% filter(Freq != 1)\n",
    "message(paste(\"eQTls colocalized with > 1 traits:\", length(unique(b$Var1))))\n",
    "summary(b$Freq)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 161,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "# total loci: 183\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. \n",
       "   1.00    4.00   12.00   26.88   35.00  227.00 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "# loci with 1 SNP: 21\n",
      "\n",
      "# loci 2-10 SNP: 63\n",
      "\n",
      "# loci >10 SNP: 99\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. \n",
       "   11.0    18.0    33.0    46.3    55.0   227.0 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "b = gwas_coloc[gwas_coloc$eqtl_id %in% a$eqtl_id,]\n",
    "message(paste(\"# total loci:\", length(unique(b$gwas_locus_id))))\n",
    "summary(b$cs_size)\n",
    "message(paste(\"# loci with 1 SNP:\", length(unique(b[b$cs_size == 1,]$gwas_locus_id))))\n",
    "message(paste(\"# loci 2-10 SNP:\", length(unique(b[b$cs_size >= 2 & b$cs_size <= 10,]$gwas_locus_id))))\n",
    "message(paste(\"# loci >10 SNP:\", length(unique(b[b$cs_size >= 10,]$gwas_locus_id))))\n",
    "\n",
    "# Mean of SNPs per CS with > 10 SNPs\n",
    "summary(b[b$cs_size >= 10,]$cs_size)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Combinatorial"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 280,
   "metadata": {},
   "outputs": [],
   "source": [
    "gwas_coloc = fread(\"reviews/tables/Table_S14_GWAS_Colocalization.txt\", data.table = F)\n",
    "eqtl_annot = fread(\"reviews/tables/Table_S11_eQTL_Annotations.txt\", data.table = F)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 164,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "# total modules: 2832\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "\n",
       "alternative_splicing      gene_expression \n",
       "                 980                 1852 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "a = eqtl_annot %>% filter(!category_annotation %like% \"singleton\" & !category_annotation %like% \"failed\" & category_annotation != \"ambiguous\")\n",
    "b = a %>% select(module_id, eqtl_phenotype, category_annotation) %>% distinct()\n",
    "message(paste(\"# total modules:\", length(unique(b$module_id))))\n",
    "table(b$eqtl_phenotype)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 173,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "# modules with GWAS: 89 3.143 2832\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "89"
      ],
      "text/latex": [
       "89"
      ],
      "text/markdown": [
       "89"
      ],
      "text/plain": [
       "[1] 89"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "\n",
       "alternative_splicing      gene_expression \n",
       "                  32                   57 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Total GWAS loci colocalized with module: 129\n",
      "\n"
     ]
    }
   ],
   "source": [
    "b = gwas_coloc[gwas_coloc$module_id %in% a$module_id,] %>% select(module_id, eqtl_phenotype, category_annotation) %>% distinct()\n",
    "message(paste(\"# modules with GWAS:\", length(unique(b$module_id)), signif(length(unique(b$module_id)) / length(unique(a$module_id)) * 100, 4), length(unique(a$module_id))))\n",
    "length(unique(b$module_id))\n",
    "table(b$eqtl_phenotype)\n",
    "\n",
    "b = gwas_coloc[gwas_coloc$module_id %in% a$module_id,]\n",
    "message(paste(\"Total GWAS loci colocalized with module:\", length(unique(b$gwas_locus_id))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 179,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "89"
      ],
      "text/latex": [
       "89"
      ],
      "text/markdown": [
       "89"
      ],
      "text/plain": [
       "[1] 89"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "# ppc-unique modules: 5\n",
      "\n",
      "# fetal-adult modules: 36\n",
      "\n",
      "# adult-only modules: 48\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "\n",
       "         adult-shared           fetal-adult           fetal-islet \n",
       "                   21                    21                     4 \n",
       " fetal-whole-pancreas       ipsc_ppc-unique          islet-unique \n",
       "                   11                     5                    22 \n",
       "whole pancreas-unique \n",
       "                    5 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "b = gwas_coloc[gwas_coloc$module_id %in% a$module_id,] %>% select(module_id, eqtl_phenotype, category_annotation) %>% distinct()\n",
    "length(unique(b$module_id))\n",
    "message(paste(\"# ppc-unique modules:\", nrow(b[b$category_annotation == \"ipsc_ppc-unique\",])))\n",
    "message(paste(\"# fetal-adult modules:\", nrow(b[b$category_annotation %in% c(\"fetal-adult\", \"fetal-islet\", \"fetal-whole-pancreas\"),])))\n",
    "message(paste(\"# adult-only modules:\", nrow(b[b$category_annotation %in% c(\"adult-shared\", \"islet-unique\", \"whole pancreas-unique\"),])))\n",
    "table(b$category_annotation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 181,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "alternative_splicing \n",
       "                   5 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "table(b[b$category_annotation == \"ipsc_ppc-unique\",]$eqtl_phenotype)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 190,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "# GWAS loci: 129\n",
      "# modules: 89\n",
      "\n",
      "# total module: 89\n",
      "\n",
      "# ppc eqtls / genes: 49 41\n",
      "\n",
      "# islet eqtls / genes: 98 75\n",
      "\n",
      "# pancreas eqtls/ genes: 71 69\n",
      "\n"
     ]
    }
   ],
   "source": [
    "b = gwas_coloc[gwas_coloc$module_id %in% a$module_id,]\n",
    "message(paste(\"# GWAS loci:\", length(unique(b$gwas_locus_id))), appendLF = F)\n",
    "message(paste(\"# modules:\", length(unique(b[b$eqtl_type %like% \"comb\",]$module_id))))\n",
    "\n",
    "c = eqtl_annot %>% filter(module_id %in% b$module_id)\n",
    "\n",
    "message(paste(\"# total module:\", length(unique(c$module_id))))\n",
    "message(paste(\"# ppc eqtls / genes:\", length(unique(c[c$tissue == \"ipsc-ppc\",]$eqtl_id)), length(unique(c[c$tissue == \"ipsc-ppc\",]$gene_id))))\n",
    "message(paste(\"# islet eqtls / genes:\", length(unique(c[c$tissue == \"islet\",]$eqtl_id)), length(unique(c[c$tissue == \"islet\",]$gene_id))))\n",
    "message(paste(\"# pancreas eqtls/ genes:\", length(unique(c[c$tissue == \"pancreas\",]$eqtl_id)), length(unique(c[c$tissue == \"pancreas\",]$gene_id))))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 193,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "# GWAS loci: 129\n",
      "\n",
      "# GWAS loci with 1 SNP: 15\n",
      "\n",
      "# GWAS loci with 2-10 SNP: 54\n",
      "\n",
      "# GWAS loci with >10 SNP: 60\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. \n",
       "  11.00   16.75   22.50   31.77   42.00  106.00 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "b = gwas_coloc[gwas_coloc$module_id %in% a$module_id & gwas_coloc$used_to_finemap == T,]\n",
    "message(paste(\"# GWAS loci:\", length(unique(b$gwas_locus_id))))\n",
    "message(paste(\"# GWAS loci with 1 SNP:\", length(unique(b[b$cs_size == 1,]$gwas_locus_id))))\n",
    "message(paste(\"# GWAS loci with 2-10 SNP:\", length(unique(b[b$cs_size >= 2 & b$cs_size <= 10,]$gwas_locus_id))))\n",
    "message(paste(\"# GWAS loci with >10 SNP:\", length(unique(b[b$cs_size > 10,]$gwas_locus_id))))\n",
    "\n",
    "# Avg snps for CS with > 10 snps\n",
    "summary(b[b$cs_size > 10,]$cs_size)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Spatiotemporally informed eQTL resource provides novel insights into GWAS signals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 297,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "   fetal-adult    fetal-islet fetal-pancreas \n",
       "            21              4             11 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Total: 36\n",
      "\n",
      "Same: 30 83.3333333333333\n",
      "\n",
      "Partial: 1 2.77777777777778\n",
      "\n",
      "Different: 5 13.8888888888889\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "\n",
       " A  B  C  D  E \n",
       "11 19  1  4  1 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "a = gwas_coloc %>% filter(eqtl_type == \"combinatorial\" & category_annotation != \"ambiguous\" & !category_annotation %like% \"fail\")\n",
    "a = a %>% filter(category_annotation %in% c(\"fetal-islet\", \"fetal-whole-pancreas\", \"fetal-adult\"))\n",
    "\n",
    "mod1 = fread(\"reviews/tables/Table_S12_ModuleSummary_Gene.txt\", data.table = F) \n",
    "mod2 = fread(\"reviews/tables/Table_S12_ModuleSummary_AltSplicing.txt\", data.table = F)\n",
    "mod = rbind(mod1, mod2)\n",
    "\n",
    "b1 = mod %>% filter(module_id %in% a$module_id)\n",
    "table(b1$category_annotation)\n",
    "\n",
    "message(paste(\"Total:\", length(unique(b1$module_id))))\n",
    "message(paste(\"Same:\", length(unique(b1[b1$egene_overlap_category %in% c(\"A\", \"B\"),]$module_id)), length(unique(b1[b1$egene_overlap_category %in% c(\"A\", \"B\"),]$module_id)) / length(unique(b$module_id)) * 100))\n",
    "message(paste(\"Partial:\", length(unique(b1[b1$egene_overlap_category %in% c(\"C\"),]$module_id)), length(unique(b1[b1$egene_overlap_category %in% c(\"C\"),]$module_id)) / length(unique(b1$module_id)) * 100))\n",
    "message(paste(\"Different:\", length(unique(b1[b1$egene_overlap_category %in% c(\"D\", \"E\"),]$module_id)), length(unique(b1[b1$egene_overlap_category %in% c(\"D\", \"E\"),]$module_id)) / length(unique(b1$module_id)) * 100))\n",
    "\n",
    "table(b1$egene_overlap_category)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# eASQTLs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 215,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "# easQTLs per tissue:\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "\n",
       "ipsc-ppc    islet pancreas \n",
       "    3959     4939     2077 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "eqtl_annot = fread(\"reviews/tables/Table_S11_eQTL_Annotations.txt\", data.table = F) \n",
    "eqtl_annot_as = eqtl_annot %>% filter(eqtl_phenotype %like% \"alt\")\n",
    "\n",
    "message(\"# easQTLs per tissue:\")\n",
    "table(eqtl_annot_as$tissue)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 222,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "# easQTL pairs\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "4868"
      ],
      "text/latex": [
       "4868"
      ],
      "text/markdown": [
       "4868"
      ],
      "text/plain": [
       "[1] 4868"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "pairs = fread(\"reviews/tables/Table_S10_GraphInput_AltSplicing.txt\", data.table = F)\n",
    "message(\"# easQTL pairs\")\n",
    "nrow(pairs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 213,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "# easQTL singletons\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table class=\"dataframe\">\n",
       "<caption>A data.frame: 3 × 2</caption>\n",
       "<thead>\n",
       "\t<tr><th scope=col>Var1</th><th scope=col>Freq</th></tr>\n",
       "\t<tr><th scope=col>&lt;fct&gt;</th><th scope=col>&lt;int&gt;</th></tr>\n",
       "</thead>\n",
       "<tbody>\n",
       "\t<tr><td>ipsc_ppc singleton</td><td> 631</td></tr>\n",
       "\t<tr><td>islet singleton   </td><td>1522</td></tr>\n",
       "\t<tr><td>pancreas singleton</td><td> 431</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "A data.frame: 3 × 2\n",
       "\\begin{tabular}{ll}\n",
       " Var1 & Freq\\\\\n",
       " <fct> & <int>\\\\\n",
       "\\hline\n",
       "\t ipsc\\_ppc singleton &  631\\\\\n",
       "\t islet singleton    & 1522\\\\\n",
       "\t pancreas singleton &  431\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "A data.frame: 3 × 2\n",
       "\n",
       "| Var1 &lt;fct&gt; | Freq &lt;int&gt; |\n",
       "|---|---|\n",
       "| ipsc_ppc singleton |  631 |\n",
       "| islet singleton    | 1522 |\n",
       "| pancreas singleton |  431 |\n",
       "\n"
      ],
      "text/plain": [
       "  Var1               Freq\n",
       "1 ipsc_ppc singleton  631\n",
       "2 islet singleton    1522\n",
       "3 pancreas singleton  431"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "message(\"# easQTL singletons\")\n",
    "data.frame(table(eqtl_annot_as[eqtl_annot_as$category_annotation %like% \"single\",]$category_annotation))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 231,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "# easQTL modules\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "980"
      ],
      "text/latex": [
       "980"
      ],
      "text/markdown": [
       "980"
      ],
      "text/plain": [
       "[1] 980"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "980"
      ],
      "text/latex": [
       "980"
      ],
      "text/markdown": [
       "980"
      ],
      "text/plain": [
       "[1] 980"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "980"
      ],
      "text/latex": [
       "980"
      ],
      "text/markdown": [
       "980"
      ],
      "text/plain": [
       "[1] 980"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "eqtl_annot_as_pass = eqtl_annot_as %>% filter(!category_annotation %like% \"ambi\" & !category_annotation %like% \"failed\")\n",
    "\n",
    "message(\"# easQTL modules\")\n",
    "length(unique(eqtl_annot_as_pass[!eqtl_annot_as_pass$category_annotation %like% \"single\",]$module_id))\n",
    "length(unique(eqtl_annot_as_pass[!eqtl_annot_as_pass$eqtl_type %like% \"single\",]$module_id))\n",
    "\n",
    "mod = fread(\"reviews/tables/Table_S12_ModuleSummary_AltSplicing.txt\", data.table = F)  %>% filter(!category_annotation %like% \"module_failed\" & !category_annotation %like% \"ambi\")\n",
    "nrow(mod)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 236,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "# tissue-unique easQTL modules\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "344"
      ],
      "text/latex": [
       "344"
      ],
      "text/markdown": [
       "344"
      ],
      "text/plain": [
       "[1] 344"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "\n",
       "ipsc_ppc-unique    islet-unique pancreas-unique \n",
       "            124             203              17 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "# easQTL in tissue-unique easQTL modules\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "\n",
       "ipsc-ppc    islet pancreas \n",
       "     266      452       37 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "message(\"# tissue-unique easQTL modules\")\n",
    "length(unique(eqtl_annot_as_pass[eqtl_annot_as_pass$category_annotation %like% \"unique\",]$module_id))\n",
    "table(mod[mod$category_annotation %like% \"unique\",]$category_annotation)\n",
    "\n",
    "message(\"# easQTL in tissue-unique easQTL modules\")\n",
    "a = eqtl_annot_as_pass[eqtl_annot_as_pass$category_annotation %like% \"unique\",]\n",
    "table(a$tissue)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 241,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "# adult-shared easQTL modules\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "225"
      ],
      "text/latex": [
       "225"
      ],
      "text/markdown": [
       "225"
      ],
      "text/plain": [
       "[1] 225"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "225"
      ],
      "text/latex": [
       "225"
      ],
      "text/markdown": [
       "225"
      ],
      "text/plain": [
       "[1] 225"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "message(\"# adult-shared easQTL modules\")\n",
    "length(unique(eqtl_annot_as_pass[eqtl_annot_as_pass$category_annotation == \"adult-shared\",]$module_id))\n",
    "length(unique(mod[mod$category_annotation == \"adult-shared\",]$module_id))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 260,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "# modules shared between iPSC-PPC and adult\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "\n",
       "   fetal-adult    fetal-islet fetal-pancreas \n",
       "           214            139             58 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "411"
      ],
      "text/latex": [
       "411"
      ],
      "text/markdown": [
       "411"
      ],
      "text/plain": [
       "[1] 411"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "\n",
       "   fetal-adult    fetal-islet fetal-pancreas \n",
       "           214            139             58 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "411"
      ],
      "text/latex": [
       "411"
      ],
      "text/markdown": [
       "411"
      ],
      "text/plain": [
       "[1] 411"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "# easQTLs per tissue in the 411 modules\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "\n",
       "ipsc-ppc    islet pancreas \n",
       "     802      561      318 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "message(\"# modules shared between iPSC-PPC and adult\")\n",
    "a = eqtl_annot_as_pass[eqtl_annot_as_pass$category_annotation %in% c(\"fetal-islet\", \"fetal-pancreas\", \"fetal-adult\"),] %>% select(category_annotation, module_id) %>% distinct()\n",
    "table(a$category_annotation)\n",
    "length(unique(a$module_id))\n",
    "\n",
    "\n",
    "a = mod[mod$category_annotation %in% c(\"fetal-islet\", \"fetal-pancreas\", \"fetal-adult\"),]\n",
    "table(a$category_annotation)\n",
    "nrow(a)\n",
    "\n",
    "message(\"# easQTLs per tissue in the 411 modules\")\n",
    "a = eqtl_annot_as_pass[eqtl_annot_as_pass$category_annotation %in% c(\"fetal-islet\", \"fetal-pancreas\", \"fetal-adult\"),]\n",
    "table(a$tissue)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 267,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "# iPSC-PPC-unique easQTLs\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "897"
      ],
      "text/latex": [
       "897"
      ],
      "text/markdown": [
       "897"
      ],
      "text/plain": [
       "[1] 897"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "\n",
       "combinatorial     singleton \n",
       "          266           631 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "# iPSC-PPC-unique easQTLs shared with adult\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "802"
      ],
      "text/latex": [
       "802"
      ],
      "text/markdown": [
       "802"
      ],
      "text/plain": [
       "[1] 802"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "# iPSC-PPC failed criteria\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "2260"
      ],
      "text/latex": [
       "2260"
      ],
      "text/markdown": [
       "2260"
      ],
      "text/plain": [
       "[1] 2260"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "message(\"# iPSC-PPC-unique easQTLs\")\n",
    "a = eqtl_annot_as_pass %>% filter(category_annotation %like% \"single\" | category_annotation %like% \"unique\") %>% filter(tissue == \"ipsc-ppc\")\n",
    "nrow(a)\n",
    "\n",
    "table(a$eqtl_type)\n",
    "\n",
    "message(\"# iPSC-PPC-unique easQTLs shared with adult\")\n",
    "a = eqtl_annot_as_pass %>% filter(!category_annotation %like% \"single\"  & !category_annotation %like% \"unique\") %>% filter(tissue == \"ipsc-ppc\")\n",
    "nrow(a)\n",
    "\n",
    "message(\"# iPSC-PPC failed criteria\")\n",
    "a = eqtl_annot_as %>% filter(category_annotation %like% \"ambi\" | category_annotation %like% \"failed\") %>% filter(tissue == \"ipsc-ppc\")\n",
    "nrow(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Characterization of fetal-adult-shared eASQTL modules in iPSC-PPC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 279,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "   fetal-adult    fetal-islet fetal-pancreas \n",
       "           214            139             58 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "411"
      ],
      "text/latex": [
       "411"
      ],
      "text/markdown": [
       "411"
      ],
      "text/plain": [
       "[1] 411"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "\n",
       "  A   B   C   D   E \n",
       "149  85  93  57  27 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "177"
      ],
      "text/latex": [
       "177"
      ],
      "text/markdown": [
       "177"
      ],
      "text/plain": [
       "[1] 177"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "\n",
       "ipsc-ppc    islet pancreas \n",
       "     386      361      184 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "a = mod[mod$category_annotation %in% c(\"fetal-islet\", \"fetal-pancreas\", \"fetal-adult\"),]\n",
    "table(a$category_annotation)\n",
    "nrow(a)\n",
    "table(a$egene_overlap_category)\n",
    "\n",
    "nrow(a[a$egene_overlap_category %in% c(\"C\", \"D\", \"E\"),])\n",
    "\n",
    "a = eqtl_annot_as_pass %>% filter(module_id %in% a[a$egene_overlap_category %in% c(\"C\", \"D\", \"E\"),]$module_id)\n",
    "table(a$tissue)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Methods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 307,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "# total gwas loci that colocalized with eQTL\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "312"
      ],
      "text/latex": [
       "312"
      ],
      "text/markdown": [
       "312"
      ],
      "text/plain": [
       "[1] 312"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "# loci by eqtl type\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "\n",
       "combinatorial     singleton \n",
       "          129           183 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "gwas_coloc = fread(\"reviews/tables/Table_S14_GWAS_Colocalization.txt\", data.table = F) %>% filter(!category_annotation %like% \"ambiguous\" & !category_annotation %like% \"failed\")\n",
    "\n",
    "message(\"# total gwas loci that colocalized with eQTL\")\n",
    "length(unique(gwas_coloc$gwas_locus_id))\n",
    "\n",
    "message(\"# loci by eqtl type\")\n",
    "a = gwas_coloc %>% select(gwas_locus_id, eqtl_type) %>% distinct()\n",
    "table(a$eqtl_type)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 315,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "# gwas loci in cred sets\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "312"
      ],
      "text/latex": [
       "312"
      ],
      "text/markdown": [
       "312"
      ],
      "text/plain": [
       "[1] 312"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "check if there are non-overlapping loci\n",
      "\n",
      "Warning message in cbind(parts$left, chars$ellip_h, parts$right, deparse.level = 0L):\n",
      "“number of rows of result is not a multiple of vector length (arg 2)”\n",
      "Warning message in cbind(parts$left, chars$ellip_h, parts$right, deparse.level = 0L):\n",
      "“number of rows of result is not a multiple of vector length (arg 2)”\n",
      "Warning message in cbind(parts$left, chars$ellip_h, parts$right, deparse.level = 0L):\n",
      "“number of rows of result is not a multiple of vector length (arg 2)”\n",
      "Warning message in cbind(parts$left, chars$ellip_h, parts$right, deparse.level = 0L):\n",
      "“number of rows of result is not a multiple of vector length (arg 2)”\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table class=\"dataframe\">\n",
       "<caption>A data.frame: 0 × 21</caption>\n",
       "<thead>\n",
       "\t<tr><th scope=col>gwas_locus_id</th><th scope=col>eqtl_id</th><th scope=col>trait_id</th><th scope=col>description</th><th scope=col>snp</th><th scope=col>pvalues.df1</th><th scope=col>MAF.df1</th><th scope=col>N.df1</th><th scope=col>V.df1</th><th scope=col>z.df1</th><th scope=col>⋯</th><th scope=col>lABF.df1</th><th scope=col>pvalues.df2</th><th scope=col>MAF.df2</th><th scope=col>N.df2</th><th scope=col>V.df2</th><th scope=col>z.df2</th><th scope=col>r.df2</th><th scope=col>lABF.df2</th><th scope=col>internal.sum.lABF</th><th scope=col>SNP.PP.H4</th></tr>\n",
       "\t<tr><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>⋯</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;int&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th></tr>\n",
       "</thead>\n",
       "<tbody>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "A data.frame: 0 × 21\n",
       "\\begin{tabular}{lllllllllllllllllllll}\n",
       " gwas\\_locus\\_id & eqtl\\_id & trait\\_id & description & snp & pvalues.df1 & MAF.df1 & N.df1 & V.df1 & z.df1 & ⋯ & lABF.df1 & pvalues.df2 & MAF.df2 & N.df2 & V.df2 & z.df2 & r.df2 & lABF.df2 & internal.sum.lABF & SNP.PP.H4\\\\\n",
       " <chr> & <chr> & <chr> & <chr> & <chr> & <dbl> & <dbl> & <int> & <dbl> & <dbl> & ⋯ & <dbl> & <dbl> & <dbl> & <int> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl> & <dbl>\\\\\n",
       "\\hline\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "A data.frame: 0 × 21\n",
       "\n",
       "| gwas_locus_id &lt;chr&gt; | eqtl_id &lt;chr&gt; | trait_id &lt;chr&gt; | description &lt;chr&gt; | snp &lt;chr&gt; | pvalues.df1 &lt;dbl&gt; | MAF.df1 &lt;dbl&gt; | N.df1 &lt;int&gt; | V.df1 &lt;dbl&gt; | z.df1 &lt;dbl&gt; | ⋯ ⋯ | lABF.df1 &lt;dbl&gt; | pvalues.df2 &lt;dbl&gt; | MAF.df2 &lt;dbl&gt; | N.df2 &lt;int&gt; | V.df2 &lt;dbl&gt; | z.df2 &lt;dbl&gt; | r.df2 &lt;dbl&gt; | lABF.df2 &lt;dbl&gt; | internal.sum.lABF &lt;dbl&gt; | SNP.PP.H4 &lt;dbl&gt; |\n",
       "|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|\n",
       "\n"
      ],
      "text/plain": [
       "     gwas_locus_id eqtl_id trait_id description snp pvalues.df1 MAF.df1 N.df1\n",
       "     V.df1 z.df1 ⋯ lABF.df1 pvalues.df2 MAF.df2 N.df2 V.df2 z.df2 r.df2\n",
       "     lABF.df2 internal.sum.lABF SNP.PP.H4"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning message in cbind(parts$left, chars$ellip_h, parts$right, deparse.level = 0L):\n",
      "“number of rows of result is not a multiple of vector length (arg 2)”\n",
      "Warning message in cbind(parts$left, chars$ellip_h, parts$right, deparse.level = 0L):\n",
      "“number of rows of result is not a multiple of vector length (arg 2)”\n",
      "Warning message in cbind(parts$left, chars$ellip_h, parts$right, deparse.level = 0L):\n",
      "“number of rows of result is not a multiple of vector length (arg 2)”\n",
      "Warning message in cbind(parts$left, chars$ellip_h, parts$right, deparse.level = 0L):\n",
      "“number of rows of result is not a multiple of vector length (arg 2)”\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<table class=\"dataframe\">\n",
       "<caption>A data.frame: 0 × 29</caption>\n",
       "<thead>\n",
       "\t<tr><th scope=col>gwas_locus_id</th><th scope=col>eqtl_id</th><th scope=col>transcript_id</th><th scope=col>gene_id</th><th scope=col>gene_name</th><th scope=col>eqtl_phenotype</th><th scope=col>tissue</th><th scope=col>eqtl_type</th><th scope=col>module_id</th><th scope=col>trait_id</th><th scope=col>⋯</th><th scope=col>likely_model</th><th scope=col>topsnp</th><th scope=col>topsnp_pp</th><th scope=col>topsnp_gwas_pval</th><th scope=col>topsnp_eqtl_pval</th><th scope=col>islet_egene_overlap</th><th scope=col>pancreas_egene_overlap</th><th scope=col>category_annotation</th><th scope=col>used_to_finemap</th><th scope=col>cs_size</th></tr>\n",
       "\t<tr><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>⋯</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;dbl&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;chr&gt;</th><th scope=col>&lt;lgl&gt;</th><th scope=col>&lt;int&gt;</th></tr>\n",
       "</thead>\n",
       "<tbody>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "A data.frame: 0 × 29\n",
       "\\begin{tabular}{lllllllllllllllllllll}\n",
       " gwas\\_locus\\_id & eqtl\\_id & transcript\\_id & gene\\_id & gene\\_name & eqtl\\_phenotype & tissue & eqtl\\_type & module\\_id & trait\\_id & ⋯ & likely\\_model & topsnp & topsnp\\_pp & topsnp\\_gwas\\_pval & topsnp\\_eqtl\\_pval & islet\\_egene\\_overlap & pancreas\\_egene\\_overlap & category\\_annotation & used\\_to\\_finemap & cs\\_size\\\\\n",
       " <chr> & <chr> & <chr> & <chr> & <chr> & <chr> & <chr> & <chr> & <chr> & <chr> & ⋯ & <chr> & <chr> & <dbl> & <dbl> & <dbl> & <chr> & <chr> & <chr> & <lgl> & <int>\\\\\n",
       "\\hline\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "A data.frame: 0 × 29\n",
       "\n",
       "| gwas_locus_id &lt;chr&gt; | eqtl_id &lt;chr&gt; | transcript_id &lt;chr&gt; | gene_id &lt;chr&gt; | gene_name &lt;chr&gt; | eqtl_phenotype &lt;chr&gt; | tissue &lt;chr&gt; | eqtl_type &lt;chr&gt; | module_id &lt;chr&gt; | trait_id &lt;chr&gt; | ⋯ ⋯ | likely_model &lt;chr&gt; | topsnp &lt;chr&gt; | topsnp_pp &lt;dbl&gt; | topsnp_gwas_pval &lt;dbl&gt; | topsnp_eqtl_pval &lt;dbl&gt; | islet_egene_overlap &lt;chr&gt; | pancreas_egene_overlap &lt;chr&gt; | category_annotation &lt;chr&gt; | used_to_finemap &lt;lgl&gt; | cs_size &lt;int&gt; |\n",
       "|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|\n",
       "\n"
      ],
      "text/plain": [
       "     gwas_locus_id eqtl_id transcript_id gene_id gene_name eqtl_phenotype\n",
       "     tissue eqtl_type module_id trait_id ⋯ likely_model topsnp topsnp_pp\n",
       "     topsnp_gwas_pval topsnp_eqtl_pval islet_egene_overlap\n",
       "     pancreas_egene_overlap category_annotation used_to_finemap cs_size"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "message(\"# gwas loci in cred sets\")\n",
    "cred_sets = fread(\"reviews/tables/Table_S15_99CredibleSet_GWAS.txt\", data.table = F)\n",
    "length(unique(cred_sets$gwas_locus_id))\n",
    "\n",
    "gwas_coloc = fread(\"reviews/tables/Table_S14_GWAS_Colocalization.txt\", data.table = F) %>% filter(!category_annotation %like% \"ambiguous\" & !category_annotation %like% \"failed\")\n",
    "\n",
    "message(\"check if there are non-overlapping loci\")\n",
    "cred_sets %>% filter(!gwas_locus_id %in% gwas_coloc$gwas_locus_id)\n",
    "gwas_coloc %>% filter(!gwas_locus_id %in% cred_sets$gwas_locus_id)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R 4.2.1",
   "language": "R",
   "name": "ir421"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.2.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
